Lattice model

This notebook simulates the evolution of the system with local dispersal, featuring three variants of the rock-paper-scissors game. All lattice-based model were developed using the Mesa framework. With Mesa's visualization tools, we can observe the spatial structures that emerge from the interactions between the agents.

Outline

• 1. Frean and Abraham model
  • 1.1. Implementation details
  • 1.2. Simulations
• 2. Model with simultaneous activation of agents
  • 2.1. Implementation details
  • 2.2. Simulations
• 3. Reichenbach, Mobilia and Frey model
  • 3.1. Implementation details
  • 3.2. Simulations
• 4. References

Importing the necessary libraries:

import sys
import os
import solara
import mesa
import numpy as np
import pandas as pd
from matplotlib.figure import Figure
import plotly.express as px
import plotly.figure_factory as ff
import matplotlib.pyplot as plt
from mesa.experimental import JupyterViz
sys.path.append(os.path.abspath('..'))
from lattice_models.model_rand_act import RSPRandAct
from lattice_models.model_sim_act import RSPSimAct
from lattice_models.model_mobility import RSPMobility
%matplotlib inline

/Users/irenetesta/anaconda3/envs/stochpy/lib/python3.11/site-packages/solara/validate_hooks.py:122: UserWarning: /Users/irenetesta/anaconda3/envs/stochpy/lib/python3.11/site-packages/mesa/experimental/jupyter_viz.py:163: JupyterViz: `use_state` found within a nested function created on line 155
To suppress this check, replace the line with:
        grid_layout, set_grid_layout = solara.use_state(grid_layout_initial)  # noqa: SH104

Make sure you understand the consequences of this, by reading about the rules of hooks at:
    https://solara.dev/documentation/advanced/understanding/rules-of-hooks

  warnings.warn(str(e))

1. Frean and Abraham model

In this model the N sites are taken to be sites in a periodic, square lattice and each agent is activated once per time-step, in random order. During each activation, an agent competes with a randomly chosen neighboring agent. If the neighboring agent can be defeated, the agent wins the competition with a probability determined by the invasion rate. If the agent wins the competition, the neighboring agent is transformed into an individual of the same species of the winning agent.

For more details, see the original paper [1].

1.1. Implementation details

The model is defined by the class RSPRandAct in the file ../lattice_models/model_rand_act.py:

%pycat ../lattice_models/model_rand_act.py

import mesa
from .patch_rand_act import PatchRandAct
from mesa.datacollection import DataCollector

class RSPRandAct(mesa.Model):
    """
    A system with three species in a competitive loop: a rock beats a pair of scissors,
    scissors beat a sheet of paper and paper beats a rock.
    """

    # key: species, value: list of species that the key species can beat
    # e.g. ROCK (0) beats SCISSOR (1), SCISSOR (1) beats PAPER (2), PAPER (2) beats ROCK (0)
    rules3 = {0: [1], 1: [2], 2: [0]}
    rules4 = {0: [1], 1: [2], 2: [3], 3: [0]}
    rules5 = {0: [1,2], 1: [2,3], 2: [3,4], 3: [4,0], 4: [0,1]}

    def __init__(self,
                 init0, init1, init2, init3, init4,
                 invrate0, invrate1, invrate2, invrate3, invrate4,
                 n_species,
                 color_map,
                 increase_rate=False,
                 width=50, height=50):
        """
        Create a new playing area of (width, height) patches.
        """
        super().__init__()

        self.n_species = n_species
        self.color_map = color_map
        self.increase_rate = increase_rate
        if self.n_species == 3:
            self.init_probabilities = [init0, init1, init2]
            self.invasion_rates = [invrate0, invrate1, invrate2]
            self.rules = self.rules3
        elif self.n_species == 4:
            self.init_probabilities = [init0, init1, init2, init3]
            self.invasion_rates = [invrate0, invrate1, invrate2, invrate3]
            self.rules = self.rules4
        else: # n_species == 5
            self.init_probabilities = [init0, init1, init2, init3, init4]
            self.invasion_rates = [invrate0, invrate1, invrate2, invrate3, invrate4]
            self.rules = self.rules5

        # set up agent scheduler
        self.schedule = mesa.time.RandomActivation(self)

        # use a simple grid, where edges wrap around.
        self.grid = mesa.space.SingleGrid(width, height, torus=True)

        # place a patch at each location, randomly initializing it according to the given initial proportions
        for _, (x, y) in self.grid.coord_iter():
            patch_init_state = self.random.choices(range(self.n_species), weights=self.init_probabilities, k=1)[0]
            patch = PatchRandAct(pos=(x, y), model=self, init_state=patch_init_state)
            self.grid.place_agent(patch, (x, y))
            self.schedule.add(patch)

        self.running = True
        
        # collect data
        model_reporter = {}
        for i in range(self.n_species):
            model_reporter[i] = lambda model, species=i: model.count_patches(species)
        if self.increase_rate:
           model_reporter['increased_rate'] = lambda model: model.invasion_rates[0] # TODO: lambda model: model.compute_average_invrate0()=
        self.datacollector = DataCollector(
            model_reporter
        )
        self.datacollector.collect(self)

    def count_patches(self, species):
        """
        Helper method to count the number of patches in the given state.
        """
        return sum([1 for patch in self.schedule.agents if patch.state == species])
    
    def compute_average_invrate0(self): # TODO: delete?
        rates = [patch.invasion_rate for patch in self.schedule.agents if patch.state == 0]
        return sum(rates) / len(rates)

    def step(self):
        """
        Have the scheduler advance each patch by one step.
        """
        self.schedule.step()
        self.datacollector.collect(self)

        n_existint_species = sum([1 for i in range(self.n_species) if self.count_patches(i) == 0])
        if n_existint_species == self.n_species-1:
            self.running = False

• The model can be simulated with 3, 4, or 5 species. When there are 4 species, species <1,3> and <2,4> do not compete with each other.
• Each agent is activated once per step, in random order, with the order reshuffled each step, through the mesa.time.RandomActivation scheduler.
• Agents are situated on a rectangular, periodic grid (i.e. torus), where each cell contains exactly at most one agent. Such a grid is created by the mesa.space.SingleGrid class.
• The grid is randomly initialized by placing agents in each cell according to the initial proportions of the species (the parameters init* are used as weights in the random.choice function).

The behavior of each agent (grid patch) is defined by the class PatchRandAct in the file ../lattice_models/patch_rand_act.py:

%pycat ../lattice_models/patch_rand_act.py

import mesa

class PatchRandAct(mesa.Agent):
    """Represents a single patch in the simulation."""

    def __init__(self, pos, model, init_state):
        """
        Create a patch, in the given state, at the given x, y position.
        """
        super().__init__(pos, model)
        self.rules = model.rules
        self.n_species = model.n_species
        self.color_map = model.color_map
        self.x, self.y = pos
        self.state = init_state
        self.invasion_rate = self.model.invasion_rates[self.state] # TODO: delete?
        self._nextState = None

    def step(self):
        """
        Compute the next state of the patch according to the rules of the game
        and the invasion rates of the species.
        """
        # get the neighbors
        neighbors = self.model.grid.get_neighbors((self.x, self.y), moore=True)

        # select a random neighbor
        neighbor = self.random.choice(neighbors)
        
        # if the neighbor can be defeated
        if neighbor.state in self.rules[self.state]:
            # change the neighbor state according to the invasion rate
            win_rate = self.model.invasion_rates[self.state]
            new_neigh_state = self.random.choices(
                [self.state, neighbor.state],
                weights=[win_rate, 1-win_rate],
                k=1)[0]
            neighbor.state = new_neigh_state
            
            # increase species 0 invasion rate if needed
            if (self.model.increase_rate and self.state==0 and \
                new_neigh_state==self.state and self.model.schedule.steps>100):
                rand = self.random.uniform(0, 1e-4)
                new_rate = self.model.invasion_rates[self.state] + rand
                if 0 < new_rate < 1:
                    self.model.invasion_rates[self.state] = new_rate

    def alternative_step(self):
        """
        Compute the next state of the patch according to the rules of the game
        and the invasion rates of the species.
        """
        # get the neighbors
        neighbors = self.model.grid.get_neighbors((self.x, self.y), moore=True)

        # select a random neighbor
        neighbor = self.random.choice(neighbors)
        
        # if the neighbor can be defeated
        if neighbor.state in self.rules[self.state]:
            # change the neighbor state according to the invasion rate
            win_rate = self.invasion_rate
            new_neigh_state = self.random.choices(
                [self.state, neighbor.state],
                weights=[win_rate, 1-win_rate],
                k=1)[0]
            if new_neigh_state != neighbor.state:
                if self.state==0:
                    neighbor.invasion_rate = self.invasion_rate
                else: 
                    neighbor.invasion_rate = self.model.invasion_rates[neighbor.state]
            neighbor.state = new_neigh_state
            
            # increase species 0 individuals invasion rate if needed
            if (self.model.increase_rate and self.state==0 and \
                new_neigh_state==self.state and self.model.schedule.steps>100):
                rand = self.random.uniform(-1e-2, 1e-2)
                new_rate = self.invasion_rate + rand
                if 0 < new_rate < 1:
                    self.invasion_rate = new_rate

• Each patch has a Moore neighborhood, i.e. the eight cells surrounding it.
• At each step, a neighbor is randomly chosen and the agent in the patch interacts with it. According to the invasion rates, the neighboring agent may be defeated and may become the same species as the agent in the patch.
• If the parameter increase_rate is set to True, after 100 steps the invasion rate of species 0 is increased whenever it replicates onto a new site. The increment is a random number uniformly distributed between 0 and 1e-4 and the new invasion rate is accepted if it is in the range $0<P_r<1$).

1.2. Simulations

Defining functions to run simulations and visualize the evolution of the system:

def run_model(model_class, params, steps):
    model = model_class(**params)
    for i in range(steps):
        model.step()
    return model

def agent_portrayal(cell):
    '''
    This function defines the appearance of the agents in the visualization.
    '''
    return {"color": cell.color_map[cell.state], "marker": "s", "size": 1}

def draw_grid(grid, agent_portrayal, title=None, ax=None):
    '''
    This function draws the lattice grid.
    '''
    def portray(g):
        x = []
        y = []
        s = []  # size
        c = []  # color
        for i in range(g.width):
            for j in range(g.height):
                content = g._grid[i][j]
                if not content:
                    continue
                if not hasattr(content, "__iter__"):
                    # Is a single grid
                    content = [content]
                for agent in content:
                    data = agent_portrayal(agent)
                    x.append(i)
                    y.append(j)
                    if "size" in data:
                        s.append(data["size"])
                    if "color" in data:
                        c.append(data["color"])
        out = {"x": x, "y": y}
        # This is the default value for the marker size, which auto-scales
        # according to the grid area.
        out["s"] = (180 / max(g.width, g.height)) ** 2
        if len(s) > 0:
            out["s"] = s
        if len(c) > 0:
            out["c"] = c
        return out

    if ax is None:
        space_fig = Figure()
        space_ax = space_fig.subplots()
    else:
        space_ax = ax
    space_ax.set_xlim(-1, grid.width)
    space_ax.set_ylim(-1, grid.height)
    space_ax.scatter(**portray(grid))
    space_ax.set_title(title)
    space_ax.axis('off')
    if ax is None:
        return space_fig

def line_plot(model, title=None, ax=None):
    '''
    This function plots the fraction of individuals of each type over time.
    '''
    #all_colors = []
    all_species = []
    patches_types = model.n_species

    if hasattr(model, 'expected_swap'):
        #all_colors.append('black')
        all_species.append('empty')
        patches_types += 1

    #all_colors.extend(['tab:red', 'purple', 'gold', 'tab:blue', 'tab:green'])
    all_species.extend(['$n_r$', '$n_s$', '$n_p$', '$n_u$', '$n_v$'])

    if ax is None:
        fig = Figure()
        axs = fig.subplots()
    else:
        axs = ax
    model_data = model.datacollector.get_model_vars_dataframe()
    
    tot_individuals = 0
    for i in range(patches_types):
        tot_individuals += model_data[i]
    for i in range(patches_types):
        model_data[i] = model_data[i] / tot_individuals
    
    model_data[[i for i in range(patches_types)]].plot(
        ax=axs,
        style='-o',
        color=[model.color_map[i] for i in range(patches_types)]
    )
    axs.legend([all_species[i] for i in range(patches_types)])
    axs.set_xlabel('Time Step')
    axs.set_ylabel('Fraction of Individuals')
    axs.set_title(title)
    
    if ax is None:
        solara.FigureMatplotlib(fig)
        return fig

def plot_rate_species0(model, title=None):
    '''
    This function plots the invasion rate of species 0 over time as well as the fraction of individuals of species 0.
    '''
    fig = Figure()
    ax = fig.subplots()
    model_data = model.datacollector.get_model_vars_dataframe()
    tot_individuals = model_data[0] + model_data[1] + model_data[2]
    model_data[0] = model_data[0] / tot_individuals
    model_data[[0, 'increased_rate']].plot(ax=ax, style='-o', color=['tab:red', 'black'])
    ax.legend(['$n_r$', 'average $P_r$'])
    ax.set_xlabel('Time Step')
    ax.set_title(title)
    solara.FigureMatplotlib(fig)
    return fig

def plot_explored_area(model, title=None):
    '''
    This function plots the explored area over time.
    '''
    fig = Figure()
    ax = fig.subplots()
    model_data = model.datacollector.get_model_vars_dataframe()
    model_data['mean_explored_area'].plot(ax=ax, style='-o', color='black')
    ax.set_xlabel('Time Step')
    ax.set_ylabel('Mean Explored Area')
    ax.set_title(title)
    solara.FigureMatplotlib(fig)
    return fig

def ternary_plot_species_evolution(df, title=None, show_markers=True):
    '''
    Plot the evolution of species proportions over time in a ternary plot.
    
    Parameters
    ----------
    df : pd.DataFrame
        A DataFrame containing the species proportions over time.
    title : str
        The title of the plot.
    '''
    fig = px.scatter_ternary(
        df,
        a=df.columns[2],
        b=df.columns[0],
        c=df.columns[1],
        color=df.index,
        color_continuous_scale='blues',
        size_max=10,
        title=title
    )
    fig.update_layout(
        width=750,
    )
    if show_markers:
        fig.update_traces(mode='lines+markers', line=dict(color='black', width=1),
            marker=dict(symbol='circle', line=dict(width=1, color='black')))
    else:
        fig.update_traces(mode='lines', line=dict(color='black', width=1))
    fig.update_layout(coloraxis_colorbar=dict(title='Time'))
    fig.show('png')

Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites and invasion rates $P_r=0.2, P_s=0.5, P_p=0.3$. The grid will be initialized with the sites being randomly assigned to one of the three species ($n_r=0.33, n_s=0.33, n_p=0.33$). In all the simulations, species R will be represented in red, species S in purple and species P in yellow.

color_map = {
    0: 'red',
    1: 'purple',
    2: 'gold',
}
model_params1 = {
    "height": 200,
    "width": 200,
    "n_species": 3,
    "init0": 0.33,
    "init1": 0.33,
    "init2": 0.33,
    "init3": 0,
    "init4": 0,
    "invrate0": 0.2,
    "invrate1": 0.5,
    "invrate2": 0.3,
    "invrate3": 0,
    "invrate4": 0,
    "color_map": color_map
}

Running the simulation for 300 steps:

model1 = run_model(model_class=RSPRandAct, params=model_params1, steps=300)
results_df = model1.datacollector.get_model_vars_dataframe().rename(columns={0: '$n_r$', 1: '$n_s$', 2: '$n_p$'})

Defining a function to build the title for the plots using the model's parameters:

def build_title(params, steps=None):
    '''
    This function builds the title for the plots using the parameters of the model.
    '''
    title = rf"$N={params['height']} \times {params['width']},"

    if 'swap_exp' in params:
        title += f" empty={params['init0']:.2f}, n_r^0={params['init1']:.2f}, n_s^0={params['init2']:.2f}, n_p^0={params['init3']:.2f}"
        if params['init4'] > 0:
            title += f", n_u^0={params['init4']:.2f}"
        if params['init5'] > 0:
            title += f", n_v^0={params['init5']:.2f}"
    else:
        title += f" n_r^0={params['init0']:.2f}, n_s^0={params['init1']:.2f}, n_p^0={params['init2']:.2f}"
        if params['init3'] > 0:
            title += f", n_u^0={params['init3']:.2f}"
        if params['init4'] > 0:
            title += f", n_v^0={params['init4']:.2f}"

    if 'invrate0' in params:
        title += f", P_r={params['invrate0']}, P_s={params['invrate1']}, P_p={params['invrate2']}"
    if 'invrate3' in params and params['invrate3'] > 0:
        title += f", P_u={params['invrate3']}"
    if 'invrate4' in params and params['invrate4'] > 0:
        title += f", P_v={params['invrate4']}"

    if 'swap_exp' in params:
        title += f", swap\_exp={params['swap_exp']}, reproduce\_exp={params['reproduce_exp']}, select\_exp={params['select_exp']}"

    if steps:
        title += f", step={steps}$"
    else:
        title += "$"
    
    return title

Plotting the evolution of the system:

ternary_plot_species_evolution(results_df, title=build_title(params=model_params1))

Unlike the long-range dispersal model, where the weakest species outcompetes the others in the same setting, here the dynamics stabilize, and the population densities converge towards a fixed point.

Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites and equal invasion rates ($Pr_=0.33, P_s=0.33, P_p=0.33$). The grid will be initialized with the sites being randomly assigned to each of the three species in fixed-point proportions ($n_r=0.33, n_s=0.33, n_p=0.33$).

model_params2 = {
    "height": 200,
    "width": 200,
    "n_species": 3,
    "init0": 0.33,
    "init1": 0.33,
    "init2": 0.33,
    "init3": 0,
    "init4": 0,
    "invrate0": 0.33,
    "invrate1": 0.33,
    "invrate2": 0.33,
    "invrate3": 0,
    "invrate4": 0,
    "color_map": color_map
}

Creating the widget to interactively simulate the system:

page = JupyterViz(
    RSPRandAct,
    model_params2,
    measures=[line_plot],
    agent_portrayal=agent_portrayal,
)
page

Cannot show widget. You probably want to rerun the code cell above (Click in the code cell, and press Shift+Enter ⇧+↩).

Since the widget cannot be saved in the notebook, let's run a simulation on a square lattice with $N=200\times 200$ sites and plot the grid after 100 steps and the evolution of the species proportions over time.

model2 = run_model(model_class=RSPRandAct, params=model_params2, steps=100)

draw_grid(
    model2.grid,
    agent_portrayal,
    title=build_title(params=model_params2, steps=100)
)

line_plot(
    model2,
    title=build_title(params=model_params2)
)

With equal invasion rates the populations form clumps.

Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites and invasion rates $P_r=0.1, P_s=0.1, P_p=0.8$. The grid is initialized with the sites being randomly assigned to each of the three species in fixed-point proportions ($n_r=0.1, n_s=0.8, n_p=0.1$).

model_params3 = {
    "height": 200,
    "width": 200,
    "color_map": color_map,
    "n_species": 3,
    "init0": 0.1,
    "init1": 0.8,
    "init2": 0.1,
    "init3": 0,
    "init4": 0,
    "invrate0": 0.1,
    "invrate1": 0.1,
    "invrate2": 0.8,
    "invrate3": 0,
    "invrate4": 0,
}

Creating the widget to interactively simulate the system:

page = JupyterViz(
    RSPRandAct,
    model_params3,
    measures=[line_plot],
    agent_portrayal=agent_portrayal,
)
page

Cannot show widget. You probably want to rerun the code cell above (Click in the code cell, and press Shift+Enter ⇧+↩).

Running a simulation and plotting the grid after 100 steps and the evolution of the species proportions over time.

model3 = run_model(model_class=RSPRandAct, params=model_params3, steps=100)

draw_grid(
    model3.grid,
    agent_portrayal,
    title=build_title(params=model_params3, steps=100)
)

line_plot(
    model3,
    title=build_title(params=model_params3)
)

When species S invasion rates is much larger than the other two, species R and P form disconnected islands amongst species S.

Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites and invasion rates $P_r=0.05, P_s=0.475, P_p=0.475$. The grid is initialized with the sites being randomly assigned to each of the three species in fixed-point proportions ($n_r=0.475, n_s=0.475, n_p=0.05$).

model_params4 = {
    "height": 200,
    "width": 200,
    "color_map": color_map,
    "n_species": 3,
    "init0": 0.475,
    "init1": 0.475,
    "init2": 0.05,
    "init3": 0,
    "init4": 0,
    "invrate0": 0.05,
    "invrate1": 0.475,
    "invrate2": 0.475,
    "invrate3": 0,
    "invrate4": 0,
}

Creating the widget to interactively simulate the system:

page = JupyterViz(
    RSPRandAct,
    model_params4,
    measures=[line_plot],
    agent_portrayal=agent_portrayal,
)
page

Cannot show widget. You probably want to rerun the code cell above (Click in the code cell, and press Shift+Enter ⇧+↩).

Running a simulation and plotting the grid after 100 steps and the evolution of the species proportions over time.

model4_100 = run_model(model_class=RSPRandAct, params=model_params4, steps=100)

draw_grid(
    model4_100.grid,
    agent_portrayal,
    title=build_title(params=model_params4, steps=100)
    )

line_plot(
    model4_100,
    title=build_title(params=model_params4)
)

Running a simulation and plotting the grid after 300 steps and the evolution of the species proportions over time.

model4_300 = run_model(model_class=RSPRandAct, params=model_params4, steps=300)

draw_grid(
    model4_300.grid,
    agent_portrayal,
    title=build_title(params=model_params4, steps=300)
    )

line_plot(
    model4_300,
    title=build_title(params=model_params4)
)

If species R invasion rate is much smaller than the other two, then species R and S have similar patchy distributions, with species P persisting in fast-moving thin fronts.

Defining the parameters to perform simulations for a range of invasion probabilities chosen to sum to unity, initializing the population densities with equal proportions:

# defining a grid of points for invasion probabilities
R, S = np.mgrid[0:1:10j, 0:1:10j]
R, S = R.ravel(), S.ravel()
P = 1 - R - S
I = np.array([R, S, P]).T
# keeping only invasion probabilities greater than 0.01
I = I[np.all(I > 0.01, axis=1)]
R = I[:, 0]
S = I[:, 1]
P = I[:, 2]
# parameters for the batch run
params_batch = {
    "height": 100,
    "width": 100,
    "color_map": color_map,
    "n_species": 3,
    "init0": 0.33,
    "init1": 0.33,
    "init2": 0.33,
    "init3": 0,
    "init4": 0,
    "invrate3": 0,
    "invrate4": 0,
}

Setting a macro to actually run the simulations or load the results of th already executed simulations from a file (executing the simulations may take a long time):

EXECUTE = False

Performing simulations:

if EXECUTE:
    # dataframe to store the results
    batch_results_df = pd.DataFrame()

    # iterating over the invasion probabilities
    for i, (Pr, Ps, Pp) in enumerate(I):
        params_batch['invrate0'] = Pr
        params_batch['invrate1'] = Ps
        params_batch['invrate2'] = Pp
        print(f"[{i+1}/{len(I)}] Running model with parameters: {params_batch}")
        # simulating with the current invasion probabilities
        results = mesa.batch_run(
            RSPRandAct,
            params_batch,
            iterations=1,
            max_steps=300,
            number_processes=1,
            data_collection_period=1,
            display_progress=True,
        )
        # concatenating the results of the current simulation
        curr_results_df = pd.DataFrame(results)
        curr_results_df['simulation_id'] = i
        tot_individuals = curr_results_df[0] + curr_results_df[1] + curr_results_df[2]
        curr_results_df.rename(columns={0: 'individuals_r', 1: 'individuals_s', 2: 'individuals_p'}, inplace=True)
        curr_results_df['$n_r$'] = curr_results_df['individuals_r'] / tot_individuals
        curr_results_df['$n_s$'] = curr_results_df['individuals_s'] / tot_individuals
        curr_results_df['$n_p$'] = curr_results_df['individuals_p'] / tot_individuals
        batch_results_df = pd.concat([batch_results_df, curr_results_df])

    # saving the results to a csv file
    batch_results_df.to_csv("../results/lattice_simulations_results.csv", index=False)
else:
    batch_results_df = pd.read_csv("../results/lattice_simulations_results.csv")

Plotting a countour summarizing the results of 28 simulations that are each 300 epochs long on a square lattice with $N=100\cdot 100$ sites and showing the mean value of $n_s$ over the final 100 epochs of each model run. Data are only shown for the region that had all its invasion probabilities greater than 0.01 as outside this region, with such a small invasion probability, the grid is too small to stabilize fluctuations in the population densities.

ns_final = np.array(batch_results_df[batch_results_df['Step']>200].groupby('simulation_id')['$n_s$'].mean())
fig = ff.create_ternary_contour(
    np.array([I.T[2], I.T[1], I.T[0]]),
    ns_final,
    pole_labels=['$P_p$', '$P_s$', '$P_r'],
    interp_mode='ilr',
    ncontours=15,
    colorscale='Portland',
    coloring=None,
    showmarkers=True,
    showscale=True,
    title=r'$\text{Mean value of }n_s\text{ over the final 100 steps of each simulation}\\\text{(population densities initialized with equal proportions)}$',
)
fig.update_layout(width=600)
fig.show('png')

The value of $n_s$ at the mean-field fixed point varies linearly from the base of the triangle to the top. Deviations from the mean-field approximation are seen toward the corners of the triangle, but $n_s$ still rises with decreasing $P_s$.

Despite the large clumps that may develop when the invasion rates are unequal, the population densities remain close to the fixed point of the mean-field equations over much of the parameter space of the invasion probabilities and the paradox, that the most aggressive species does not have the largest population, still holds.

Let's define the parameters to perform a simulation of a system with 3 species in which the invasion rate of a species (species R) is increased over time (increase_rate=True), on a square lattice with $N=200\cdot 200$ sites and initial invasion rates $P_r=0.05, P_s=0.5, P_p=0.3$. The grid is initialized with the sites being randomly assigned to each of the three species in fixed-point proportions ($n_r=0.5, n_s=0.3, n_p=0.05$).

model_params5 = {
    "height": 200,
    "width": 200,
    "color_map": color_map,
    "n_species": 3,
    "init0": 0.5,
    "init1": 0.3,
    "init2": 0.05,
    "init3": 0,
    "init4": 0,
    "invrate0": 0.05,
    "invrate1": 0.5,
    "invrate2": 0.3,
    "invrate3": 0,
    "invrate4": 0,
    "increase_rate": True
}

Creating the widget to interactively simulate the system:

page = JupyterViz(
    RSPRandAct,
    model_params5,
    measures=[line_plot, plot_rate_species0],
    agent_portrayal=agent_portrayal,
)
page

Cannot show widget. You probably want to rerun the code cell above (Click in the code cell, and press Shift+Enter ⇧+↩).

model5 = run_model(model_class=RSPRandAct, params=model_params5, steps=400)

line_plot(
    model5,
    title=build_title(params=model_params5)
)

plot_rate_species0(
    model5,
    title=build_title(params=model_params5)
)

An increase in the invasion rate of species R results in a decrease of its population.

2. Model with simultaneous activation of agents

This is a variant of the model where the N sites are again taken to be sites in a periodic, square lattice, but at each time step, all agents are activated simultaneously. During activation, each agent competes with all its neighboring agents. If an agent is surrounded by at least three neighbors of the species that defeat it, it is transformed into an individual of that species. Note that invasion rates are not considered in this model.

2.1. Implementation details

The model is defined by the class RSPSimAct in the file ../lattice_models/model_sim_act.py:

%pycat ../lattice_models/model_sim_act.py

import mesa
from .patch_sim_act import PatchSimAct
from mesa.datacollection import DataCollector

class RSPSimAct(mesa.Model):
    """
    A system with three species in a competitive loop: a rock beats a pair of scissors,
    scissors beat a sheet of paper and paper beats a rock.
    """

    rules3 = {0: [1], 1: [2], 2: [0]}
    rules4 = {0: [1], 1: [2], 2: [3], 3: [0]}
    rules5 = {0: [1,2], 1: [2,3], 2: [3,4], 3: [4,0], 4: [0,1]}

    def __init__(self, init0, init1, init2, init3, init4, n_species, color_map, width=50, height=50):
        """
        Create a new playing area of (width, height) patches.
        """
        super().__init__()

        self.n_species = n_species
        self.color_map = color_map

        if self.n_species == 3:
            self.probabilities = [init0, init1, init2]
            self.threshold = 3
            self.rules = self.rules3
        elif self.n_species == 4:
            self.probabilities = [init0, init1, init2, init3]
            self.threshold = 2
            self.rules = self.rules4
        else: # n_species == 5
            self.probabilities = [init0, init1, init2, init3, init4]
            self.threshold = 3
            self.rules = self.rules5

        # set up agent scheduler
        self.schedule = mesa.time.SimultaneousActivation(self)

        # use a simple grid, where edges wrap around.
        self.grid = mesa.space.SingleGrid(width, height, torus=True)

        # place a patch at each location, initializing it as ROCK, SCISSOR, or PAPER
        for _, (x, y) in self.grid.coord_iter():
            patch_init_state = self.random.choices(range(self.n_species), weights=self.probabilities, k=1)[0]
            patch = PatchSimAct(pos=(x, y), model=self, init_state=patch_init_state)
            self.grid.place_agent(patch, (x, y))
            self.schedule.add(patch)

        self.running = True
        
        # collect data
        model_reporter = {}
        for i in range(self.n_species):
            model_reporter[i] = lambda model, species=i: model.count_patches(species)
        self.datacollector = DataCollector(
            model_reporter
        )
        self.datacollector.collect(self)

    def count_patches(self, species):
        """
        Helper method to count the number of patches in the given state.
        """
        return sum([1 for patch in self.schedule.agents if patch.state == species])

    def step(self):
        """
        Have the scheduler advance each patch by one step.
        """
        self.schedule.step()
        self.datacollector.collect(self)

        n_existint_species = sum([1 for i in range(self.n_species) if self.count_patches(i) == 0])
        if n_existint_species == self.n_species-1:
            self.running = False

• The model can be simulated with 3, 4, or 5 species. When there are 4 species, species <1,3> and <2,4> do not compete with each other.
• Simultaneous activation of agents is implemented through the mesa.time.SimultaneousActivation scheduler.
• Agents are situated on a rectangular, periodic grid (i.e. torus), where each cell contains exactly at most one agent. Such a grid is created by the mesa.space.SingleGrid class.
• The grid is randomly initialized by placing agents in each cell according to the initial proportions of the species (the parameters init* are used as weights in the random.choice function).

The behavior of each agent (grid patch) is defined by the class PatchSimAct in the file ../lattice_models/patch_sim_act.py:

%pycat ../lattice_models/patch_sim_act.py

import mesa
import numpy as np

class PatchSimAct(mesa.Agent):
    """Represents a single patch in the simulation."""

    def __init__(self, pos, model, init_state):
        """
        Create a patch, in the given state, at the given x, y position.
        """
        super().__init__(pos, model)
        self.rules = model.rules
        self.n_species = model.n_species
        self.color_map = model.color_map
        self.threshold = model.threshold
        self.x, self.y = pos
        self.state = init_state
        self._nextState = None

    def step(self):
        """
        Compute the next state of the path according to the rules of the game and
        the number of defeats of each type.
        The state is not changed here, but is just computed and stored in self._nextState,
        because our current state may still be necessary for our neighbors
        to calculate their next state.
        """

        # assume nextState is unchanged, unless changed below
        self._nextState = self.state

        # get the neighbors
        neighbors = self.model.grid.get_neighbors((self.x, self.y), moore=True)

        # count the number of defeats of each type
        defeat_counts = np.zeros(self.n_species)
        for neighbor in neighbors:
            if neighbor.state != self.state and self.state in self.rules[neighbor.state]:
                defeat_counts[neighbor.state] += 1
        
        # compute the winning species
        winners = np.argwhere(defeat_counts == np.max(defeat_counts)).reshape(-1)
        self.random.shuffle(winners)
        for i in range(len(winners)):
            if defeat_counts[winners[i]] >= self.threshold:
                self._nextState = winners[i]
                break

    def advance(self):
        """
        Set the state to the new state, computed in the method step().
        """
        self.state = self._nextState

• Each patch has a Moore neighborhood, i.e. the eight cells surrounding it.
• At each step, an agent interacts with all its neighbors. If an agent is surrounded by at least three neighbors of the species that defeat it, it is transformed into an individual of that species. In case of a tie (possible with 5 species), the new species is chosen randomly among the species that defeated the agent.

2.2. Simulations

Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites, initialized with the sites being randomly assigned to each of the three species in equal proportions (init0=0.33, init1=0.33, init2=0.33).

model_params6 = {
    "height": 200,
    "width": 200,
    "n_species": 3,
    "init0": 0.33,
    "init1": 0.33,
    "init2": 0.33,
    "init3": 0,
    "init4": 0,
    "color_map": color_map
}

Creating the widget to interactively simulate the system:

page = JupyterViz(
    RSPSimAct,
    model_params6,
    measures=[line_plot],
    agent_portrayal=agent_portrayal,
)
page

Cannot show widget. You probably want to rerun the code cell above (Click in the code cell, and press Shift+Enter ⇧+↩).

Running a simulation and plotting the grid after 200 steps and the evolution of the species proportions over time.

model6 = run_model(model_class=RSPSimAct, params=model_params6, steps=200)

draw_grid(
    model6.grid,
    agent_portrayal,
    title=build_title(params=model_params6, steps=200)
)

line_plot(
    model6,
    title=build_title(params=model_params6)
)

results_df = model6.datacollector.get_model_vars_dataframe().rename(columns={0: '$n_r$', 1: '$n_s$', 2: '$n_p$'})
ternary_plot_species_evolution(results_df, title=build_title(params=model_params6))

The system dynamics reach a fixed point where the population densities are equal and it leads to the emergence of entangled rotating spiral waves.

Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites, initialized with the sites being randomly assigned to each of the three species according to the following probabilities init0=0.2, init1=0.4, init2=0.4.

model_params7 = {
    "height": 200,
    "width": 200,
    "n_species": 3,
    "init0": 0.2,
    "init1": 0.4,
    "init2": 0.4,
    "init3": 0,
    "init4": 0,
    "color_map": color_map
}

Creating the widget to interactively simulate the system:

page = JupyterViz(
    RSPSimAct,
    model_params7,
    measures=[line_plot],
    agent_portrayal=agent_portrayal,
)
page

Cannot show widget. You probably want to rerun the code cell above (Click in the code cell, and press Shift+Enter ⇧+↩).

Running a simulation and plotting the grid after 400 steps and the evolution of the species proportions over time.

model7 = run_model(model_class=RSPSimAct, params=model_params7, steps=400)

draw_grid(
    model7.grid,
    agent_portrayal,
    title=build_title(params=model_params7, steps=400)
)

line_plot(
    model7,
    title=build_title(params=model_params7)
)

results_df = model7.datacollector.get_model_vars_dataframe().rename(columns={0: '$n_r$', 1: '$n_s$', 2: '$n_p$'})
ternary_plot_species_evolution(results_df, title=build_title(params=model_params7))

The system eventually stabilizes at a fixed point with equal population densities, resulting in the formation of intricate octagonal patterns.

Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites, initialized with the sites being randomly assigned to each of the three species according to the following probabilities init0=0.1, init1=0.3, init2=0.5.

model_params8 = {
    "height": 200,
    "width": 200,
    "n_species": 3,
    "init0": 0.1,
    "init1": 0.3,
    "init2": 0.5,
    "init3": 0,
    "init4": 0,
    "color_map": color_map
}

Creating the widget to interactively simulate the system:

page = JupyterViz(
    RSPSimAct,
    model_params8,
    measures=[line_plot],
    agent_portrayal=agent_portrayal,
)
page

Cannot show widget. You probably want to rerun the code cell above (Click in the code cell, and press Shift+Enter ⇧+↩).

Running a simulation and plotting the evolution of the species proportions during 100 steps.

model8 = run_model(model_class=RSPSimAct, params=model_params8, steps=100)

line_plot(
    model8,
    title=build_title(params=model_params8)
)

This time the species with lower initial population density (species R) dominates the system and the other two species go extinct.

Let's define the parameters to perform a simulation of a system with 4 species, on a square lattice with $N=200\cdot 200$ sites, initialized with the sites being randomly assigned to each of the three species in equal proportions (init0=0.25, init1=0.25, init2=0.25, init3=0.25, init4=0).

model_params9 = {
    "height": 200,
    "width": 200,
    "n_species": 4,
    "init0": 0.25,
    "init1": 0.25,
    "init2": 0.25,
    "init3": 0.25,
    "init4": 0,
    "color_map": {
        0: 'red',
        1: 'purple',
        2: 'gold',
        3: 'green'
    }
}

Creating the widget to interactively simulate the system:

page = JupyterViz(
    RSPSimAct,
    model_params9,
    measures=[line_plot],
    agent_portrayal=agent_portrayal,
)
page

Cannot show widget. You probably want to rerun the code cell above (Click in the code cell, and press Shift+Enter ⇧+↩).

Running a simulation and plotting the evolution of the species proportions during 200 steps.

model9 = run_model(model_class=RSPSimAct, params=model_params9, steps=200)

draw_grid(
    model9.grid,
    agent_portrayal,
    title=build_title(params=model_params9, steps=200)
)

line_plot(
    model9,
    title=build_title(params=model_params9)
)

With 4 species, the system dynamics stabilize and regular spiral patterns emerge.

Let's define the parameters to perform a simulation of a system with 5 species, on a square lattice with $N=200\cdot 200$ sites, initialized with the sites being randomly assigned to each of the three species in equal proportions (init0=0.2, init1=0.2, init2=0.2, init3=0.2, init4=0.2).

model_params10 = {
    "height": 200,
    "width": 200,
    "n_species": 5,
    "init0": 0.2,
    "init1": 0.2,
    "init2": 0.2,
    "init3": 0.2,
    "init4": 0.2,
    "color_map": {
        0: 'red',
        1: 'purple',
        2: 'gold',
        3: 'green',
        4: 'blue',
    }
}

Creating the widget to interactively simulate the system:

page = JupyterViz(
    RSPSimAct,
    model_params10,
    measures=[line_plot],
    agent_portrayal=agent_portrayal,
)
page

Cannot show widget. You probably want to rerun the code cell above (Click in the code cell, and press Shift+Enter ⇧+↩).

Running a simulation and plotting the evolution of the species proportions during 100 steps.

model10 = run_model(model_class=RSPSimAct, params=model_params10, steps=100)

draw_grid(
    model10.grid,
    agent_portrayal,
    title=build_title(params=model_params10, steps=100)
)

line_plot(
    model10,
    title=build_title(params=model_params10)
)

With 5 species, initially distributed in equal proportions, the system oscillates and again spiral patterns emerge.

3. Reichenbach, Mobilia and Frey model

In this model the N sites are taken to be sites in a periodic, square lattice. Each patch can be occupied by one of the species or can be blank. Each tick, the following types of events happen at defined average rates:

• Select event: Two random neighbors compete with each other. The losing patch becomes blanks.
• Reproduce event: Two random neighbors attempt to reproduce. If one of the neighbors is blank, it acquires the color of the other. Nothing happens if neither neighbor is blank.
• Swap event: Two random neighbors swap color. This represents the organisms moving.

Therefore, this model combines the local competition and reproduction from the Frean and Abraham model with spatial migration, a ubiquitous feature of real ecosystems.

For more details, see the original paper [2].

The exact number of, for instance, swap events that occur each tick is drawn from a Poisson distribution with mean equal to $(number\_ of\_ patches) * 10 ^{swap-exponent}$. A Poisson distribution defines how many times a particular event occurs given an average rate for that event assuming that the occurrences of that event are independent. Here, the occurrences of the events are approximately independent since they're being performed by different organisms.

The events occur in a random order involving random pairs of neighbors.

See the original paper for more details: Reichenbach, Tobias, Mauro Mobilia, and Erwin Frey. "Mobility promotes and jeopardizes biodiversity in rock–paper–scissors games." Nature 448.7157 (2007): 1046-1049.

3.1. Implementation details

The model is defined by the class RSPMobility in the file ../lattice_models/model_mobility.py:

%pycat ../lattice_models/model_mobility.py

import mesa
import numpy as np
from .patch_mobility import PatchMobility
from mesa.datacollection import DataCollector

class RSPMobility(mesa.Model):
    """
    Represents the 2-dimensional array of patches in Rock-Scissors-Paper Game.
    """

    rules3 = {1: [2], 2: [3], 3: [1]}
    rules4 = {1: [2], 2: [3], 3: [4], 4: [1]}
    rules5 = {1: [2,3], 2: [3,4], 3: [4,5], 4: [5,1], 5: [1,2]}

    def __init__(
            self,
            init0, init1, init2, init3, init4, init5,
            swap_exp, reproduce_exp, select_exp,
            n_species, color_map, width=50, height=50
        ):
        """
        Create a new playing area of (width, height) patches.
        """
        super().__init__()

        self.n_species = n_species

        event_repetitions = (width * height) // 3
        self.expected_swap = (10**swap_exp)*event_repetitions
        self.expected_reproduce = (10**reproduce_exp)*event_repetitions
        self.expected_select = (10**select_exp)*event_repetitions

        self.color_map = color_map

        if self.n_species == 3:
            self.probabilities = [init0, init1, init2, init3]
            self.rules = self.rules3
        elif self.n_species == 4:
            self.probabilities = [init0, init1, init2, init3, init4]
            self.rules = self.rules4
        else: # n_species == 5
            self.probabilities = [init0, init1, init2, init3, init4, init5]
            self.rules = self.rules5

        # Set up the grid and schedule.

        self.schedule = mesa.time.RandomActivation(self)

        # Use a simple grid, where edges wrap around.
        self.grid = mesa.space.SingleGrid(width, height, torus=True)

        # Place a patch at each location, initializing it as ROCK, SCISSOR, or PAPER
        for _, (x, y) in self.grid.coord_iter():
            patch_init_state = self.random.choices(range(0, self.n_species+1), weights=self.probabilities, k=1)[0]
            patch = PatchMobility(pos=(x, y), model=self, init_state=patch_init_state)
            self.grid.place_agent(patch, (x, y))
            self.schedule.add(patch)

        self.running = True
        
        model_reporter = {}
        for i in range(self.n_species+1):
            model_reporter[i] = lambda model, species=i: model.count_patches(species)
        
        self.datacollector = DataCollector(
            model_reporter
        )
        self.datacollector.collect(self)

    def count_patches(self, species):
        """
        Helper method to count the number of patches in the given state.
        """
        return sum([1 for patch in self.schedule.agents if patch.state == species])

    def step(self):
        """
        Compute the number of global events and have the scheduler advance each
        patch by one step.
        """

        # create a list of events with random Poisson distribution
        self.events = []
        self.events.extend(['swap'] * int(np.random.poisson(self.expected_swap)))
        self.events.extend(['reproduce'] * int(np.random.poisson(self.expected_reproduce)))
        self.events.extend(['select'] * int(np.random.poisson(self.expected_select)))
        
        # shuffle events
        self.random.shuffle(self.events)

        self.schedule.step()
        self.datacollector.collect(self)

        n_existint_species = sum([1 for i in range(1, self.n_species+1) if self.count_patches(i) == 0])
        if n_existint_species == self.n_species-1:
            self.running = False

• The model can be simulated with 3, 4, or 5 species. init0 represents the initial proportion of empty patches. When there are 4 species, species <1,3> and <2,4> do not compete with each other.
• Each agent is activated once per step, in random order, with the order reshuffled each step, through the mesa.time.RandomActivation scheduler.
• Agents are situated on a rectangular, periodic grid (i.e. torus), where each cell contains exactly at most one agent. Such a grid is created by the mesa.space.SingleGrid class.
• The grid is randomly initialized by placing agents in each cell according to the initial proportions of the species (the parameters init* are used as weights in the random.choice function).
• This model uses an event-based approach. It calculates how many events of each type should occur each tick and then makes agents execute those events in a random order. Note that we have to compute the number of global events rather than the number of actions that each individual patch performs since the execution of those events has to be random between the patches. That is, a single patch can't perform all of their actions in one go: suppose they end up executing 5 swaps in one tick. The swaps would be shuffling things around locally rather than allowing for an organism to travel multiple steps. Hence, this code creates a list with an entry for each event type that should occur each tick, then shuffles that list so that the events are in a random order. We then randomly activate agents who pop an event off the list, and, together with random neighboring patches, run the corresponding event.

The behavior of each agent (grid patch) is defined by the class PatchMobility in the file ../lattice_models/patch_mobility.py:

%pycat ../lattice_models/patch_mobility.py

import mesa

class PatchMobility(mesa.Agent):
    """Represents a single patch in the simulation."""

    def __init__(self, pos, model, init_state):
        """
        Create a patch, in the given state, at the given x, y position.
        """
        super().__init__(pos, model)
        self.rules = model.rules
        self.n_species = model.n_species
        self.color_map = model.color_map
        self.x, self.y = pos
        self.state = init_state
        self._nextState = None

    def step(self):
        """
        Compute the next state of the patch and/or of a randomly chosen
        neighbor, according to the event that is going to happen.
        """
        # Get the neighbors
        neighbors = self.model.grid.get_neighbors((self.x, self.y), moore=False) # 4 neighbors

        # get the action
        if len(self.model.events) == 0:
            return
        action = self.model.events.pop(0)

        # select a random neighbor
        neighbor = self.random.choice(neighbors)

        # act
        if action == 'swap':
            self._nextState = neighbor.state
            neighbor.state = self.state
            self.state = self._nextState
        elif action == 'select':
            if neighbor.state != 0 and self.state in self.rules[neighbor.state]:
                self.state = 0
            if self.state != 0 and neighbor.state in self.rules[self.state]:
                neighbor.state = 0
        else:
            if neighbor.state == 0:
                neighbor.state = self.state
            elif self.state == 0:
                self.state = neighbor.state

• Each patch has a Von Neumann neighborhood, i.e. the four cells surrounding it.
• Agent state 0 represents an empty patch.
• At each step, a neighbor is randomly chosen and the agent in the patch interacts with it, according to the event type.

Defining a function to compute the percentage of each event type at each step:

def compute_mobility_perc(params):
    total = 10**(params['swap_exp']) + 10**(params['reproduce_exp']) + 10**(params['select_exp'])
    print(f"Swap %: {100*(10**(params['swap_exp'])) / total:.2f}")
    print(f"Reproduce %: {100*(10**(params['reproduce_exp'])) / total:.2f}")
    print(f"Select %: {100*(10**(params['select_exp'])) / total:.2f}")

3.2. Simulations

Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=150\cdot 150$ sites, initialized with the sites being randomly assigned to each of the three species and empty sites in equal proportions (init0=0.25, init1=0.25, init2=0.25, init3=0.25). Black patches represent empty sites.

color_map_mobility = {
    0: 'black',
    1: 'red',
    2: 'purple',
    3: 'gold'
}
model_params_mobility = {
    "height": 150,
    "width": 150,
    "n_species": 3,
    "init0": 0.25,
    "init1": 0.25,
    "init2": 0.25,
    "init3": 0.25,
    "init4": 0,
    "init5": 0,
    "color_map": color_map_mobility
}

We will run simulations setting swap exponent to -1, 0, 0.5 and 1 and fixing reproduce and select exponents to 0.

Setting parameters and computing the percentage of each event type at each step.

model_params11 = model_params_mobility.copy()
model_params11['swap_exp'] = -1
model_params11['reproduce_exp'] = 0
model_params11['select_exp'] = 0
print("Configuration 1")
compute_mobility_perc(model_params11)
print("----------------")

model_params12 = model_params_mobility.copy()
model_params12['swap_exp'] = 0
model_params12['reproduce_exp'] = 0
model_params12['select_exp'] = 0
print("Configuration 2")
compute_mobility_perc(model_params12)
print("----------------")

model_params13 = model_params_mobility.copy()
model_params13['swap_exp'] = 0.5
model_params13['reproduce_exp'] = 0
model_params13['select_exp'] = 0
print("Configuration 3")
compute_mobility_perc(model_params13)
print("----------------")

model_params14 = model_params_mobility.copy()
model_params14['swap_exp'] = 1
model_params14['reproduce_exp'] = 0
model_params14['select_exp'] = 0
print("Configuration 4")
compute_mobility_perc(model_params14)
print("----------------")

Configuration 1
Swap %: 4.76
Reproduce %: 47.62
Select %: 47.62
----------------
Configuration 2
Swap %: 33.33
Reproduce %: 33.33
Select %: 33.33
----------------
Configuration 3
Swap %: 61.26
Reproduce %: 19.37
Select %: 19.37
----------------
Configuration 4
Swap %: 83.33
Reproduce %: 8.33
Select %: 8.33
----------------

Running simulations for 300 steps.

model11 = run_model(model_class=RSPMobility, params=model_params11, steps=300)

model12 = run_model(model_class=RSPMobility, params=model_params12, steps=300)

model13 = run_model(model_class=RSPMobility, params=model_params13, steps=300)

model14 = run_model(model_class=RSPMobility, params=model_params14, steps=300)

Visualizing the grid after 300 steps:

fig, axs = plt.subplots(2, 2, figsize=(25,15))
draw_grid(
    model11.grid,
    agent_portrayal,
    title=build_title(params=model_params11, steps=300),
    ax=axs[0][0]
)
draw_grid(
    model12.grid,
    agent_portrayal,
    title=build_title(params=model_params12, steps=300),
    ax=axs[0][1]
)
draw_grid(
    model13.grid,
    agent_portrayal,
    title=build_title(params=model_params13, steps=300),
    ax=axs[1][0]
)
draw_grid(
    model14.grid,
    agent_portrayal,
    title=build_title(params=model_params14, steps=300),
    ax=axs[1][1]
)
for ax in axs.flat:
    ax.set_aspect('equal')

As in the model with simultaneous activation of agents, species self-arrange by forming patterns of moving spirals. Increasing the swap exponent (i.e. the mobility of agents), the spiral structures grow in size.

Visualizing the evolution of the system:

fig, axs = plt.subplots(2, 2, figsize=(25,10))
line_plot(
    model11,
    title=build_title(params=model_params11),
    ax=axs[0][0]
)
line_plot(
    model12,
    title=build_title(params=model_params12),
    ax=axs[0][1]
)
line_plot(
    model13,
    title=build_title(params=model_params13),
    ax=axs[1][0]
)
line_plot(
    model14,
    title=build_title(params=model_params14),
    ax=axs[1][1]
)

Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=150\cdot 150$ sites, initialized with the sites being randomly assigned to each of the three species and empty sites in equal proportions (init0=0.25, init1=0.25, init2=0.25, init3=0.25), and setting the swap exponent to 1 and the reproduce and select exponents to -0.7.

model_params15 = {
    "height": 100,
    "width": 100,
    "n_species": 3,
    "init0": 0.25,
    "init1": 0.25,
    "init2": 0.25,
    "init3": 0.25,
    "init4": 0,
    "init5": 0,
    "swap_exp": 1,
    "reproduce_exp": -0.7,
    "select_exp": -0.7,
    "color_map": color_map_mobility
}

This corresponds to the following percentages of events at each step:

compute_mobility_perc(model_params15)

Swap %: 96.16
Reproduce %: 1.92
Select %: 1.92

Running a simulation for 2000 steps.

model15 = run_model(model_class=RSPMobility, params=model_params15, steps=2000)

Visualizing the evolution of the system:

line_plot(
    model15,
    title=build_title(params=model_params15)
)

When mobility is high, biodiversity is jeopardized and lost: the system adopts a uniform state where only one species is present, while the others have died out. Which species remains is subject to a random process, all species having equal chances to survive in the model.

Let's define the parameters to perform a simulation of a system with 4 species, on a square lattice with $N=150\cdot 150$ sites, initialized with the sites being randomly assigned to each of the five species and empty sites in equal proportions and setting swap, reproduce and select exponents to 0.

model_params16 = {
    "height": 150,
    "width": 150,
    "n_species": 4,
    "init0": 0.2,
    "init1": 0.2,
    "init2": 0.2,
    "init3": 0.2,
    "init4": 0.2,
    "init5": 0,
    "swap_exp": 0,
    "reproduce_exp": 0,
    "select_exp": 0,
    "color_map": {
        0: 'black',
        1: 'red',
        2: 'purple',
        3: 'gold',
        4: 'green'
    }
}

Creating the widget to interactively simulate the system:

page = JupyterViz(
    RSPMobility,
    model_params16,
    measures=[line_plot],
    agent_portrayal=agent_portrayal,
)
page

Cannot show widget. You probably want to rerun the code cell above (Click in the code cell, and press Shift+Enter ⇧+↩).

Running a simulation and plotting the evolution of the species proportions during 300 steps.

model16 = run_model(model_class=RSPMobility, params=model_params16, steps=300)

draw_grid(
    model16.grid,
    agent_portrayal,
    title=build_title(params=model_params16, steps=300)
)

line_plot(
    model16,
    title=build_title(params=model_params16)
)

With 4 species, the system dynamics stabilize with non competing species forming islands.

Let's define the parameters to perform a simulation of a system with 5 species, on a square lattice with $N=150\cdot 150$ sites, initialized with the sites being randomly assigned to each of the five species and empty sites in equal proportions and setting swap, reproduce and select exponents to 0.

model_params17 = {
    "height": 150,
    "width": 150,
    "n_species": 5,
    "init0": 1/6,
    "init1": 1/6,
    "init2": 1/6,
    "init3": 1/6,
    "init4": 1/6,
    "init5": 1/6,
    "swap_exp": 0,
    "reproduce_exp": 0,
    "select_exp": 0,
    "color_map": {
        0: 'black',
        1: 'red',
        2: 'purple',
        3: 'gold',
        4: 'green',
        5: 'blue',
    }
}

Creating the widget to interactively simulate the system:

page = JupyterViz(
    RSPMobility,
    model_params17,
    measures=[line_plot],
    agent_portrayal=agent_portrayal,
)
page

Cannot show widget. You probably want to rerun the code cell above (Click in the code cell, and press Shift+Enter ⇧+↩).

Running a simulation and plotting the evolution of the species proportions during 300 steps.

model17 = run_model(model_class=RSPMobility, params=model_params17, steps=300)

draw_grid(
    model17.grid,
    agent_portrayal,
    title=build_title(params=model_params17, steps=300)
)

line_plot(
    model17,
    title=build_title(params=model_params17)
)

Also with 5 species, entangled rotating spiral waves emerge.

4. References

• [1] Frean, Marcus, and Edward R. Abraham. "Rock–scissors–paper and the survival of the weakest." Proceedings of the Royal Society of London. Series B: Biological Sciences 268.1474 (2001): 1323-1327.
• [2] Reichenbach, Tobias, Mauro Mobilia, and Erwin Frey. "Mobility promotes and jeopardizes biodiversity in rock–paper–scissors games." Nature 448.7157 (2007): 1046-1049.